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Abstract 

We  construct  non-local  artificial  boundary  conditions  (ABCs)  for  the  numerical  sim¬ 
ulation  of  genuinely  time-dependent  acoustic  waves  that  propagate  from  a  compact 
source  in  an  unbounded  unobstructed  space.  The  key  property  used  for  obtaining 
the  ABCs  is  the  presenee  of  laeunae,  i.e.,  sharp  aft  fronts  of  the  waves,  in  wave-type 
solutions  in  odd-dimension  spaces.  This  property  can  be  considered  a  manifestation 
of  the  Huygens’  principle.  The  ABCs  are  obtained  directly  for  the  discrete  formula¬ 
tion  of  the  problem.  They  truncate  the  original  unbounded  domain  and  guarantee 
the  complete  transparency  of  the  new  outer  boundary  for  all  the  outgoing  waves.  A 
central  feature  of  the  proposed  ABCs  is  that  the  extent  of  their  temporal  nonlocality 
is  fixed  and  limited,  and  it  does  not  come  at  the  expense  of  simplifying  the  original 
model.  It  is  rather  a  natural  consequence  of  the  existence  of  lacunae,  which  is  a  fun¬ 
damental  property  of  the  corresponding  solutions.  The  proposed  ABCs  can  be  built 
for  any  consistent  and  stable  finite-difference  scheme.  Their  accuracy  can  always 
be  made  as  high  as  that  of  the  interior  approximation,  and  it  will  not  deteriorate 
even  when  integrating  over  long  time  intervals.  Besides,  the  ABCs  are  most  flexible 
from  the  standpoint  of  geometry  and  can  handle  irregular  boundaries  on  regular 
grids  with  no  fitting/adaptation  needed  and  no  accuracy  loss  induced.  Finally,  they 
allow  for  a  wide  range  of  model  settings.  In  particular,  not  only  one  can  analyze  the 
simplest  advective  acoustics  case  with  the  uniform  background  flow,  but  also  the 
case  when  the  waves’  source  (or  scatterer)  is  engaged  in  an  accelerated  motion. 
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1  Introduction 


Artificial  Boundary  Conditions  (ABCs)  is  a  common  name  for  a  group  of 
methods  employed  for  solving  infinite-domain  problems  on  a  computer.  ABCs 
facilitate  truncation  of  the  original  unbounded  domain,  and  provide  the  re¬ 
quired  closure  for  the  resulting  truncated  formulation.  The  literature  on  the 
subject  of  ABCs  is  broad,  and  we  refer  the  reader  to  the  review  papers  [1-3]. 

Two  opposite  trends  can  be  identihed  in  building  the  ABCs.  High  accuracy 
requirements  typically  imply  that  the  ABCs  should  be  nonlocal.  In  partic¬ 
ular,  the  exact  ABCs  are  always  nonlocal  in  multidimensional  settings.  For 
unsteady  problems,  this  also  means  nonlocality  in  time.  Clearly,  nonlocality 
of  the  ABCs  may  result  in  high  computational  costs  and  elaborate  implemen¬ 
tation  strategies.  As  such,  local  alternatives  obtained  either  independently  or 
as  a  approximations  of  nonlocal  ABCs  become  viable.  These  ABCs  are  usually 
inexpensive  and  easy  to  implement,  but  may  lack  computational  accuracy. 

In  the  current  paper  we  propose  the  ABCs  for  the  numerical  simulation  of 
non-dispersive  unsteady  acoustic  waves.  We  assume  that  there  is  a  (complex) 
phenomenon/process  conhned  to  a  bounded  region  that  manifests  itself  by  the 
radiation  of  acoustic  waves  in  the  far  held.  The  waves  subsequently  propagate 
across  the  unbounded  and  unobstructed  space,  which  is  assumed  isotropic 
and  homogeneous.  The  proposed  ABCs  are  nonlocal  and  guarantee  that  the 
accuracy  of  the  boundary  treatment  can  always  be  made  at  least  as  high  as 
that  of  the  interior  discretization.  However,  the  key  feature  of  the  proposed 
ABCs  is  that  the  extent  of  their  temporal  nonlocality  is  limited,  and  does  not 
increase  as  the  time  elapses.  The  bound  on  temporal  nonlocality  comes  as  a 
consequence  of  the  presence  of  lacunae,  i.e.,  sharp  aft  fronts  of  the  waves,  in 
the  solutions  of  the  linearized  Euler  equations.  In  general,  existence  of  the 
lacunae  is  a  fundamental  property  of  wave-type  solutions  in  odd-dimension 
spaces.  It  is  often  referred  to  in  connection  with  the  Huygens’  principle. 

The  ABCs  are  constructed  by  decomposing  the  original  problem  into  the  inte¬ 
rior  and  auxiliary  sub-problems.  The  latter  is  linear  and  homogeneous  through¬ 
out  the  entire  space,  and  is  integrated  with  the  help  of  lacunae.  In  so  doing, 
the  interior  solution  is  used  to  generate  sources  for  the  auxiliary  problem,  and 
the  auxiliary  solution  is  used  to  provide  the  missing  boundary  data,  i.e.,  the 
require  closure,  for  the  interior  problem. 

Several  other  nonlocal  ABCs’  methodologies  for  unsteady  waves  have  been 
recently  advocated  in  the  literature,  most  notably  [4-10];  we  also  mention  the 
survey  paper  [11]  and  the  bibliography  there.  Compared  to  these  techniques, 
our  approach  has  a  number  of  distinctive  characteristics.  Its  high  accuracy  and 
restricted  temporal  nonlocality  have  already  been  mentioned.  Besides,  the  pro- 
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posed  lacunae-based  ABCs  are  obtained  directly  for  the  discrete  formulation 
of  the  problem,  and  can  supplement  any  consistent  and  stable  hnite-difference 
scheme.  In  other  words,  they  bypass  the  common  stages  of  hrst  deriving  and 
then  approximating  the  continuous  ABCs,  see  [2],  The  lacunae-based  ABCs 
are  especially  designed  to  withstand  long-term  numerical  integration  with  no 
deterioration  of  accuracy.  Moreover,  neither  they  are  restricted  to  any  par¬ 
ticular  shape  of  artihcial  boundary  nor  require  grid  htting.  They  also  enable 
analysis  of  a  variety  of  problem  formulations,  e.g.,  the  case  of  sound  sources 
moving  with  acceleration,  which  is  equivalent  to  advective  acoustics  with  un¬ 
steady  background  flow.  We  note  that  the  source  motion  does  not  have  to  be 
considered  relative  to  an  ambient  medium  only.  Via  the  appropriate  Galileo 
transform,  it  can  be  considered  relative  to  a  given  mean  flow  as  well. 

The  ABCs  that  we  derive  in  the  current  paper  for  acoustics  generalize  and 
further  extend  our  previous  work  [12,13],  in  which  a  similar  methodology  was 
introduced  for  the  scalar  wave  equation.  The  current  approach  also  hts  into 
the  general  theoretical  framework  of  [14]. 


2  Lacunae  of  the  Wave  Equation 


Consider  a  three-dimensional  wave  equation,  x  =  {xi,X2,X3): 


1 

dt'^ 


dxi  dxl  dxl  j 


f{x,t),  t>0, 


t=o 


0, 


(la) 

(lb) 


with  homogeneous  initial  conditions.  For  every  {x,  t),  the  solution  p  =  p{x,  t) 
of  problem  (la),  (lb)  is  given  by  the  Kirchhoff  integral  (see,  e.g.,  [15]): 


p{x,t) 


Q<Ct 


Q 


(2) 


where  ^  =  (6,6,6),  Q  =  |a;  -  ^|  =  {xi  -  6)^  +  {x2  -  6)^  +  (^3  -  6)^ 
and  =  d6'^6'^6-  If  we  assume  that  the  right-hand  side  (RHS)  f{x,t) 
of  equation  (la)  is  compactly  supported  in  space-time  on  the  domain  Q  C 
X  [0,  -I-cxd),  then  formula  (2)  immediately  implies  that 


p{x,t) 


0  for  {x,t)  G  Pi 

itr)eQ 


a;  —  ^1  <  c(t 


r),  t  >  r}  . 


(3) 


The  region  of  space-time  dehned  by  formula  (3)  is  called  lacuna  of  the  solution 
p  =  p{x,t).  This  region  is  obtained  as  the  intersection  of  all  characteristic 
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cones  of  equation  (la)  once  the  vertex  of  the  cone  sweeps  the  support  of  the 
RHS:  supp/  C  Q.  From  the  standpoint  of  physics,  lacuna  corresponds  to 
that  part  of  space-time,  on  which  the  waves  generated  by  the  sources  f{x,t), 
supp/  C  Q,  have  already  passed,  and  the  solution  has  become  zero  again. 
The  phenomenon  of  lacunae  is  inherently  three-dimensional.  The  surface  of 
the  lacuna  represents  the  trajectory  of  aft  (trailing)  fronts  of  the  waves.  The 
existence  of  aft  fronts  in  odd-dimension  spaces  is  known  as  the  Huygens’  prin¬ 
ciple,  as  opposed  to  the  so-called  wave  diffusion  which  takes  place  in  even- 
dimension  spaces,  see,  e.g.,  [15].  Let  us  also  note  that  the  aft  fronts  and  the 
lacunae  would  still  be  present  in  the  solution  of  the  wave  equation  (la)  if  the 
homogeneous  initial  conditions  (lb)  were  replaced  by  some  inhomogeneous 
initial  conditions  with  compact  support. 

The  notion  of  lacunae  (or  lacunas)  was  hrst  introduced  and  studied  by 
Petrowsky  in  [16]  (see  also  an  account  in  [17,  Chapter  VI])  for  a  variety  of 
hyperbolic  equations  and  systems;  and  general  characterization  of  their  coef- 
hcients  was  provided  that  would  guarantee  existence  of  the  lacunae.  However, 
since  work  [16]  no  constructive  examples  of  either  equations  or  systems  have 
been  obtained  for  which  lacunae  would  be  present  in  the  solutions,  besides  the 
actual  wave  equation  (la),  as  well  as  those  equations  that  either  reduce  to  or 
are  derived  from,  the  wave  equation. 


3  The  Acoustics  System  of  Equations 


The  acoustics  (linearized  Euler’s)  system  in  its  simplest  form  governs  the 
propagation  of  sound  in  an  ambient  compressible  fluid  [18,  Chapter  VHI]: 


1  dp 

+  PO  V  ■  U  =  PoQ'voh 

ot 

du 

Po-^  +  Vp  =  bvoi- 


(4) 


In  (4),  c  is  the  speed  of  sound,  and  the  variables  u  =  u{x,t)  and  p  =  p{x,t) 
represent  small  perturbations  of  the  velocity  and  pressure,  respectively.  Adi¬ 
abatic  law  p  =  c^p  has  been  applied  that  relates  p  to  the  density  perturba¬ 
tion  p.  The  quantity  po  is  the  constant  background  density.  The  source  term 
Qvoi  =  Qvoiix,  t)  that  alters  the  balance  of  mass  in  the  system  is  known  as  vol¬ 
ume  velocity  per  unit  volume,  and  the  source  term  bvoi  =  Ko\{x,t)  that  alters 
the  balance  of  momentum  is  known  as  force  per  unit  volume  (see,  e.g.,  [19]). 


Proposition  1  Assume  that  the  initial  velocity  field  is  conservative: 
u{x,0)  =  Vipoi^)-  Then,  the  velocity  potential  exists  in  the  solutions  of 
system  (4)  for  all  t  >  0  if  and  only  if  the  forcing  bvoi  is  also  conservative: 
b^o\{x,t)  =  '\/'ip{x,t).  In  this  case,  the  potential  satisfies  the  wave  equation. 
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PROOF.  Let  us  assume  that  the  force  held  is  conservative,  and  integrate  the 
second  equation  of  system  (4)  with  respect  to  time  starting  from  t  =  0: 


PQU{x,t)  =  poU{x,Q)  -  [  Vp{x,T)dT  +  [  Kol{x, T)dT 

Jo  Jo 

=  Po'^^oix)-  /  '\/p{x,T)dT  +  /  V'?/’(a;,r)dr 
Jo  Jo 


=  V 


def 


PoPo{x)  -  f  p{x,T)dT+  [  '4){x,T)dT 
Jo  Jo 


(5) 


=  poVv?(a;,f). 


Therefore,  we  obtain: 


u{x,t)  =  Vp{x,t),  = -p(a;,t) +^(a;,t),  (6) 

and  by  differentiating  the  second  relation  in  (6)  w.r.t.  time  and  subsequently 
substituting  it  into  the  hrst  equation  of  system  (4),  we  conclude  that  the 
potential  p  =  p{x,t)  dehned  by  (5)  will  satisfy  the  wave  equation  (la)  with 
the  RHS  given  by  /(a;,t)  =  -q^oi{x,t)  + 

Conversely,  let  us  assume  that  the  velocity  held  is  conservative:  u{x,t)  = 
Vp{x,  t).  Then,  the  momentum  equation  in  system  (4)  can  be  recast  as  follows: 


dp 


=  h. 


vob 


which  means  that  the  forcing  must  be  conservative:  h^o\{x,t)  =  V'ip{x,t).  □ 


Proposition  1  also  implies  that  if  b^oi{x,t)  =  0  then  the  potential  always  exists. 
In  this  case,  the  second  equality  of  (6)  reduces  to  the  conventional  relation 
between  the  potential  and  pressure:  Po%  =  —p,  see  [18,  Chapter  VIII]. 

Applying  the  gradient  operator  to  the  wave  equation  for  the  potential,  we 
conclude  that  the  velocity  vector  also  satishes  the  wave  equation: 


1  d'^u 


—  Au 


“V^vol  + 


1  dKoi 
Pqc^  dt 


(7) 


Next,  by  diherentiating  the  hrst  equation  of  (4)  with  respect  to  time,  taking 
the  divergence  of  the  second  equation  of  (4),  and  substituting  the  result  into 
the  hrst  one,  we  arrive  at  the  wave  equation  for  p\ 


1  d'^p 
dt"^ 


—  Ap 


-V- 


hvol  +  Po 


'9gvoi 
dt  ■ 


(8) 
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The  key  consideration  of  interest  is  that  if  the  RHSs  of  system  (4)  are  com¬ 
pactly  supported  in  space  and  time  on  some  domain  Q  C  x  [0,  -|-cxd),  then 
the  RHSs  of  both  equation  (7)  and  equation  (8)  will  also  be  compactly  sup¬ 
ported  on  the  same  domain  Q.  Consequently,  solutions  of  equations  (7)  and 
(8)  will  have  lacunae  of  the  same  shape  as  prescribed  by  formula  (3).  Thus, 
we  have  arrived  at  the  following 

Proposition  2  Let  the  RHSs  of  system  (4)  be  compactly  supported  in  space 
and  time:  suppq^ofxR)  C  Q  and  supp  ^  Q,  where  Q  C  x 

[0, -|-cx)).  Let  also  u{x,t)  =  Vy){xR).  Then,  solutions  of  the  acoustics  sys¬ 
tem  (4)  with  homogeneous  initial  conditions  will  have  lacunae  of  the  same 
geometry  as  provided  by  formula  (3). 

In  other  words,  we  have  shown  that  existence  of  the  velocity  potential  is  suf¬ 
ficient  for  the  acoustics  system  (4)  to  have  lacunae  in  its  solutions;  and  that 
the  geometry  of  the  lacunae  is  determined  by  that  of  the  support  of  the  RHSs, 
as  in  the  case  of  the  wave  equation.  We  do  not  know  whether  this  condition 
is  also  necessary  for  having  the  lacunae.  However,  in  the  view  of  the  comment 
provided  in  the  end  of  Section  2,  we  will  hereafter  be  using  the  foregoing  suf- 
hcient  condition  as  the  only  reliable  indication  of  the  presence  of  lacunae  in 
the  solutions  of  system  (4). 

The  group  of  numerical  algorithms  for  acoustics  that  we  describe  hereafter 
will  be  essentially  based  on  the  presence  of  lacunae.  However,  the  Kirchhoff 
formula  (2)  will  never  be  used  as  an  actual  part  of  the  algorithm  construction, 
it  will  only  be  needed  at  the  theoretical  stage,  for  determining  the  shape  of 
the  lacunae  that  will  later  be  incorporated  into  a  purely  hnite-difference  con¬ 
text.  Previously,  the  idea  of  using  the  Huygens’  principle  for  constructing  the 
ABCs  was  promoted  by  Ting  &  Miksis  [20]  and  Givoli  &  Cohen  [21].  Both 
papers,  however,  have  suggested  to  use  numerical  quadratures  to  approximate 
the  integral  (2),  and  then  couple  it  with  the  interior  solution.  Moreover,  the 
approach  of  [20]  has  never  been  actually  implemented  in  a  practical  compu¬ 
tational  setting,  whereas  the  approach  of  [21]  required  artihcial  dissipation  to 
be  added  to  the  scheme  to  fix  the  arising  instabilities. 


4  Lacunae-Based  Integration  of  the  Acoustics  System 


We  will  be  looking  for  an  irrotational  solution  of  system  (4)  on  a  bounded 
domain  S  =  S{t)  C.  We  assume  that  this  is  a  domain  of  fixed  shape 
that  has  a  hnite  diameter  d  for  all  t  >  0.  We  also  assume  that  Vt  >  0  : 
{supp  gvoi(a;,  t)  flM^}  C  S(t)  &  {supp  fovoi(a;,t)  nM^}  C  S(t).  In  other  words, 
we  want  to  compute  the  acoustic  held  on  a  given  domain  that  also  contains 
all  the  held  sources.  While  always  remaining  hnite  in  size,  this  domain  S{t)  is 
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allowed  to  move  across  the  space  according  to  a  general  law: 


Uq  =  tto(t), 


Xo  =  iCo(t)  =  [  Uo{T)dT, 

Jo 


(9) 


where  Uq  and  Xq  are  the  velocity  vector  and  coordinates  of  a  given  point  in¬ 
side  S{t),  respectively.  The  only  limitation  that  we  impose  is  that  motion  (9) 
be  subsonic:  ma,xt\tLo(t)\  =  k  <  c.  Let  now  gvoi  and  bvoi  be  some  station¬ 
ary  acoustic  sources  (RHSs  to  system  (4)),  i.e.,  compactly  supported  func¬ 
tions  of  X  on  the  hxed  domain  S'(O)  for  all  t  >  0:  {supp  gvoi(a;,  f)  fl  C 
5'(0)  &  {suppb^oiix,t)  n  C  S'(O).  Then,  we  can  incorporate  the  transla¬ 
tional  motion  (9)  into  the  acoustics  system  (4)  as  follows,  see  [12]: 


1  dp{x,t) 
dt 


+  poV  ■  u{x,t)  =  poqvo\{x 


Xo(t),t)  =  Poq^oi{x,t), 


du{x,  t) 


-h  Vp{x,t) 


Kol{x  -  Xo{t),t) 


Kol{x,t). 


(10) 


With  no  loss  of  generality,  the  initial  conditions  for  integrating  system  (10) 
will  always  be  assumed  homogeneous.  It  is  easy  to  see  from  (10)  that  the  time- 
dependent  nature  of  the  acoustic  held  is  caused  by  the  motion  of  the  sound 
sources,  on  which  the  genuine  unsteadiness  of  the  sound  generating  mecha¬ 
nisms  may  be  “superimposed.”  The  domain  of  integration  S{t)  has  a  hxed 
size/shape  and  traces  the  motion  of  the  sources.  As  shown  in  [12,  Appendix] 
using  the  Galileo  transform,  system  (10)  is  equivalent  to  advective  acoustics 
with  stationary  sources  and  unsteady  background  how  —v^  =  —ttQ(t). 


Let  us  emphasize  that  the  setup  we  have  introduced  is  quite  general  and  also 
includes  the  case  of  the  waves’  sources  that  move  relative  to  a  given  mean  how 
(rather  than  relative  to  the  ambient  huid  only).  Indeed,  a  similar  argument 
based  on  the  Galileo  transform  would  imply  that  having  the  mean  how  Vq  and 
the  source  motion  tio  on  top  of  it  is  equivalent  to  either  stationary  sources  in 
the  mean  how  'Wo  —  fp)  or  alternatively,  sources  that  move  with  the  velocity 
— ■yo  +  'fio  through  the  ambient  medium.  The  latter  setup  is  obviously  the  same 
as  (9),  as  long  as  the  condition  max^  |  —  VQ(t)  +  tio(t)  \  =  k  <  c  is  met. 


The  case  of  primary  interest  for  our  analysis  will  be  that  of  continuously 
operating  sources  in  (10),  t  G  [0,-1-cx)).  Let  us,  however,  assume  for  the  time 
being  that  not  only  the  RHSs  in  system  (10)  are  compactly  supported  in 
space,  but  also  that  their  “lifespan”  in  time  is  hnite:  suppq^oi{x,t)  G  Q  and 
supp  hvoK^^G)  ^  Qj  where  x  [0, -|-cx))  3  Q  =  {{x,t)  x  G  S{t),tQ  <t<  ti}. 
Then,  Proposition  2  implies  that  no  later  than 


,  I  d  +  -  to){c  +  k)  _ 

t  —  02  =  to  -\ - -  =  iQ  +  Jint 

C  —  k 


(11) 


all  of  the  domain  S{t)  will  fall  into  the  lacuna  dehned  by  formula  (3)  (see 
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[12,13]  for  a  more  detailed  argument),  and  will  remain  inside  the  lacuna  con¬ 
tinuously  thereafter,  i.e.,  for  all  t  >  t2-  In  other  words,  for  the  solution  of 
system  (10),  we  will  have  p{x,t)  =  0  and  u{x,t)  =  0  when  x  G  S(t)  and 
t  >t2,  where  t2  is  dehned  by  formula  (11). 

Next,  we  realize  that  during  the  time  interval  Tint  no  wave  can  travel  in 
space  further  away  than  the  distance  cTjnt  from  the  boundary  of  the  domain 
S{to).  This  means  that  we  will  also  have  p{x,t)  =  0  and  u{x,t)  =  0  for 
dist  (a;,  S'(to))  >  cTint  and  to  <  t  <  t2-  As  such,  instead  of  the  free  unob¬ 
structed  space  outside  S{t)  we  may  consider  outer  boundaries  with  arbitrary 
(reflecting)  properties.  As  long  as  none  of  these  boundaries  is  located  closer 
than  cTint  to  S(to),  the  solution  of  (10)  inside  S(t)  is  not  going  to  feel  their 
presence  for  to  <  t  <  t2. 

In  fact,  the  foregoing  limitation  for  the  location  of  outer  boundaries  can  even 
be  relaxed  if  instead  of  requiring  that  no  wave  may  reach  an  outer  boundary 
before  t  =  t2  we  introduce  a  weaker  requirement  that  no  reflected  wave  may 
reach  S{t)  before  t  =  t2-  The  latter  consideration  easily  translates  into  the 
following  estimate  for  the  minimal  distance  between  the  domain  S{to)  and  the 
allowed  location  of  any  reflecting  boundary  (see  [12,13]  for  more  detail): 


_  c  +  k 

•^min  —  r-v  -^int‘ 


(12) 


We  note  that  both  estimates  (11)  and  (12)  are  conservative  in  the  sense  that 
they  do  not  take  into  account  the  direction  of  the  motion.  In  case  the  motion 
of  the  sources  is  characterized  by  a  specihc  or  predominant  direction,  then  the 
quantity  can  be  further  reduced  in  the  orthogonal  directions. 


Altogether  we  conclude  that  the  solution  of  system  (10)  driven  by  the  sources 
compactly  supported  on  the  domain  Q  =  {{x,t)  x  G  S(t),to  <  t  <  ti}  and 
subject  to  the  homogeneous  initial  conditions,  can  be  obtained  on  S{t)  for  all 
t  >  0  as  follows.  First,  this  solution  is  obviously  equal  to  zero  for  0  <  t  <  to- 
Next,  on  the  time  interval  to  <  t  <  t2,  see  formula  (11),  system  (10)  should 
be  integrated  on  the  auxiliary  domain  of  size 


Z  —  d  -\-  —  d  -\-  (^c  -\-  k)T{YO,  (13) 

centered  around  S(to),  which  in  any  event  is  going  to  yield  the  correct  solution 
inside  S(t).  Finally,  for  all  t  >  t2  the  solution  on  S(t)  will  be  equal  to  zero 
again  because  all  the  waves  will  have  left  the  domain  by  t  =  ^2  (lacuna). 


Let  us  now  address  the  case  of  continuously  operating  sources.  For  that  pur- 


pose,  we  introduce  a  partition  of  unity  on  the  semi-infinite  interval  t  >  0: 


Vt>0:  = 

j=0 


(14) 


where  T  >  0  and  ^  <  cr  <  1  are  two  parameters,  and  0  =  0(t)  is  a  smooth, 
even,  “hat”-type  function  with  supp0(t)  C 


T  Tl. 

^  2  ’  2h 


t  >  f. 


m  = 


1,  0  <  t  <  (a  —  |)T, 

1  —  0((tT  —  t),  {a  —  l)T  <  t  <  f , 

[0(— t),  t  <  0. 


Then,  we  dehne  compactly  supported  sources  for  j  =  0, 1,  2, . . . : 

t)  =  Qvoiix,  t)Q{t  -  aTj),  supp  gS(a;,  t)  C  Qj, 
=  Ko\{x,t)e{t  -  aTj),  supp  b^{!i{x,t)  C  Qj, 
where  according  to  (15): 

Qj  =  {(®,^)|®  e  S{t),  (aj  -  ^)t  <  t  <  (aj  +  ^)t}, 
and  consider  a  collection  of  sub-problems  driven  by  the  RHSs  (16): 


1 


+  poV  ■  =  pogi^i. 


dt 


dt 


u^^\x,  t) 


t={<^3-h)T 


=  0. 


(15) 


(16) 


(17) 


(18) 


As  formula  (14)  is  a  partition  of  unity,  we  have: 

OO 

q^o\{x,  t)  =  Y,  i)  and  Ko\{.x,  t)  =  ^  ^)> 

i=o  i=o 

and  because  of  the  linear  superposition,  the  solution  p{x,  t),  u{x,  t)  of  system 
(10)  subject  to  the  homogeneous  initial  conditions  at  t  =  0  can  be  expanded 
in  terms  of  the  solutions  to  systems  (18): 


p{x,t)  =  '^p^^\x,t)  and  u{x,t)  =  '^u^^\x,t). 

j=0  j=0 


(19) 


The  series  (19)  are  formally  inhnite.  However,  for  any  f  >  0  and  x  &  S{t)  each 
will,  in  fact,  contain  only  a  hnite  hxed  number  of  nonzero  terms.  Indeed,  due 
to  the  causality  for  a  given  t  >  0  and  all  {aj  —  \)T>t,  i.e.,  all  j  >  (|i  +  |)/o', 
we  will  have  p^^\x,  t)  =  0  and  u^^x,  t)  =  0  on  the  entire  space  Moreover, 
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multiplication  by  the  function  0,  see  (16),  that  is  only  a  function  of  time  will 
not  alter  the  conservative  nature  of  6voi-  Therefore,  Proposition  2  will  apply 
to  the  solution  of  every  problem  (18),  j  =  0, 1, 2, ....  As  such,  if  we  interpret 
the  moment  {aj  —  \)T  of  the  inception  of  source  j  as  and  the  moment 
{aj  +  1)T  of  its  cessation  as  ti  \  then  according  to  formula  (11); 


d  +  T{c  +  k) 
c  —  k 


=  +  Tint 


p^^\x,t)  =  0  and  =  0  for  x  e  S{t), 


(20) 


i.e.,  the  domain  S{t)  will  be  entirely  inside  the  lacuna  starting  from  = 
Alternatively,  this  means  that  for  any  t  >  Oandallj  <  (^-^  +  |)/o', 
where  Tint  =  as  in  formula  (20),  the  terms  p^^\x,t)  and  u^^\x,t)  in 

the  series  (19)  will  be  equal  to  zero  for  x  G  S(t).  Consequently,  for  alH  >  0 
and  X  E  S(t)  we  can  replace  expansions  (19)  with 
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p{x,t)  =  ^  p^^\x,t)  and  u{x,t)  =  ^ 


m2 


U 


ij)t 


x,t), 


(21) 


j=mi 


j=mi 


where  mi  = 
integer  part. 


(^  +  l)/^ 


m2  = 


(T  +  1)/(T  +1,  and  [  ■  ]  denotes  the 


his  implies  that  for  any  t  >  0  and  x  G  S{t)  the  number  of 
terms  m  =  m2  —  mi  +  1  in  the  sum  (21),  and  as  such,  the  number  of  non- zero 
terms  in  the  expansion  (19),  may  not  exceed  +  2.  Most  important,  this 


number  m  does  not  increase  as  the  time  t  elapses,  because  the  interval  Tint 
introduced  in  (20)  depends  only  on  the  partition  size  T  for  the  sources,  the 
geometry,  the  propagation  speed,  and  the  maximum  speed  of  motion. 


Let  us  now  assume  that  each  problem  (18)  is  integrated  individually  on  an 
appropriate  domain  of  size  Z,  see  formula  (13),  by  means  of  a  consistent  and 
stable  hnite-difference  scheme.  For  every  given  j,  system  (18)  only  needs  to  be 
integrated  from  tQ^  =  {aj  —  |)T  till  +  Tint  because  for  all  subsequent 

moments  of  time  its  solution  on  S{t)  will  be  equal  to  zero.  Consequently,  the 
following  convergence  estimates  will  hold  for  ^  X  G  S{t): 

\\p^^\x,t)  -pl^\x,t)\\  <  Kjh^, 

\\u^^\x,t)  -  ui^\x,t)\\  < 

where  a  is  the  order  of  convergence,  h  denotes  the  generic  grid  size,  and  the 
functions  p^l\x,t)  and  u^\x,t)  denote  the  discrete  solution  of  system  (18) 
for  a  given  j.  The  constant  Kj  on  the  right-hand  side  of  each  inequality  (22) 
does  not  depend  on  h,  but  may  depend  on  and  as  well  as  on  Tint. 

We  emphasize  that  the  quantity  T^t  does  not  depend  on  j.  Moreover,  it  is  nat¬ 
ural  to  assume  that  the  derivatives  of  the  functions  q^^iix,  t)  and  b!fj{x,  t)  are 
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uniformly  bounded  with  respect  to  j.  In  this  case,  there  will  be  a  j-independent 
constant  K  =  K{q^oi,  bvob^int)  such  that  Vj  =  0, 1,2, :  Kj  <  K.  Then, 
using  representations  (21)  one  can  easily  transform  the  individual  convergence 
estimates  (22)  into  the  overall  temporally  uniform  grid  convergence  estimate 
for  p  and  u  that  would  hold  for  t  >  0  and  x  G  S(t): 

\\p{x,t) -ph{x,t)\\  <mKh^, 

\\u{x,t)  —  Uh{x,t)\\  <  mKh"^. 

m2  / -x  m2  / -N 

In  formulae  (23),  ph{x,t)  =  Ph  and  Uh{x,t)  =  '^h  A 

j=mi  j=mi 

detailed  proof  of  this  result  for  the  wave  equation  can  be  found  in  [12]. 

In  practical  terms,  the  temporally  uniform  grid  convergence  guaranteed  by 
estimates  (23)  means  that  accuracy  of  the  numerical  solution  of  system  (10), 
if  computed  using  lacunae,  i.e.,  by  solving  a  set  of  systems  (18)  and  then  em¬ 
ploying  representation  (21),  will  not  deteriorate  even  when  system  (10)  is  inte¬ 
grated  over  arbitrarily  long  time  intervals.  In  other  words,  one  should  expect 
that  there  will  be  no  long-time  error  buildup.  This  is,  in  fact,  a  key  distinction 
between  the  foregoing  lacunae-based  algorithm  and  traditional  time-marching 
techniques  that  may  be  applied  to  computing  the  unsteady  acoustic  helds. 

Indeed,  the  phenomenon  of  error  accumulation  during  long  runs  is  well-known 
in  the  context  of  building  computational  methods  for  time-dependent  prob¬ 
lems.  This  issue  has  been  recognized  in  the  literature  as  an  outstanding  unre¬ 
solved  question  in  numerical  PDEs  for  many  years,  since  the  first  systematic 
convergence  studies  for  discrete  approximations  have  been  conducted  in  the 
hfties.  At  the  analysis  stage,  it  manifests  itself  by  the  growth  of  the  stability 
constants  with  time.  If,  for  example,  system  (10)  needs  to  be  integrated  over 
the  interval  [0,Tfinai],  then  the  stability  constant  of  the  scheme  will,  generally 
speaking,  depend  on  TfinaO  K  =  K{  ■  ,  Tfinai),  and  will  actually  increase  with 
the  increase  of  Tfinai-  This  is,  of  course,  the  exact  same  phenomenon  as  the  de¬ 
pendence  of  Kj  on  Tint  in  formulae  (22).  The  growth  of  the  stability  constants 
with  Tfinai  is  equivalent  to  non-uniformity  of  the  grid  convergence  in  time,  and 
all  conventional  discrete  approximations  that  can  be  and  are  used  in  modern 
numerical  methods  are  known  to  suffer  from  this  dehciency.  In  practice,  this 
implies  that  any  given  approximation  can  be  used  for  only  a  limited  interval 
of  time  if  the  acceptable  level  of  error  is  prescribed.  To  advance  further  in 
time  with  no  loss  of  accuracy  a  hner  approximation  is  needed  from  the  very 
beginning,  which  obviously  prompts  the  increase  of  the  overall  computational 
cost.  The  latter  increase  may  quickly  become  prohibitive. 

Using  the  language  of  wave  physics,  one  can,  e.g.,  attribute  the  long-term  error 
buildup  to  either  numerical  dissipation,  or  dispersion  (phase  error),  or  both. 
But  no  matter  what  its  actual  mechanism  is,  it  may  result  in  an  unacceptable 
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loss  of  accuracy  by  the  solution  within  a  hnite  period  of  time. 


The  lacunae-based  algorithm  allows  us  to  circumvent  this  difficulty  due  to  the 
temporally  uniform  convergence  (23).  Moreover,  if  system  (10)  were  to  be  in¬ 
tegrated  on  the  large  interval  [0,  Tgnai]  using  a  straightforward  time-marching 
algorithm,  it  would  have  also  required  a  large  domain  in  space,  of  the  size 
roughly  2cTfi„ai.  This  is  typically  not  feasible.  On  the  other  hand,  implemen¬ 
tation  of  the  lacunae-based  algorithm  allows  us  to  perform  the  integration  on 
the  domain  of  a  hxed  and  non-increasing  size  Z  determined  by  formula  (13). 

It  is  important  to  mention  that  smoothness  plays  a  key  role  in  the  design  of  the 
lacunae-based  algorithm.  In  particular,  the  function  0(t)  of  (15)  that  helps  us 
build  the  partition  of  unity  (14)  has  to  be  chosen  sufficiently  smooth  so  that 
the  dependence  of  the  stability  constants  Kj  on  the  properties  of  individual 
RHSs  and  see  (22),  be  not  worse  than  that  in  the  original  scheme  with 
non-partitioned  source  terms.  In  this  paper,  we  leave  out  the  detailed  analysis 
that  involves  the  quantitative  smoothness  characteristics,  and  instead  refer 
the  reader  to  our  previous  work  [12]. 

Implementation  of  the  lacunae-based  algorithm  now  needs  to  be  discussed.  In 
theory,  system  (18)  for  every  given  j  =  0,1,2,...  should  be  integrated  on  its 
own  auxiliary  region  of  size  Z,  see  (13),  centered  around  the  domain  S(tQ^), 
where  =  {aj  —  l)T  and  the  location  of  S{tg  ^)  is,  in  turn,  determined  by 
the  reference  point  a^(f[/^),  see  formula  (9).  However,  it  is  more  convenient 
to  consider  periodic  boundary  conditions  with  the  period  Z  in  all  coordinate 
directions.  In  this  case,  motion  (9)  should  be  interpreted  as  the  motion  on  a 
three-dimensional  toroidal  surface,  and  all  spatial  locations  shall  be  converted 
to  the  periodic  setting:  x  x  =  {xi,X2,xg),  where  Xi  =  Xi  —  [^]Z,  i  = 
1,  2,  3.  In  so  doing,  all  systems  (18)  can  basically  be  solved  on  one  and  the 
same  domain  with  periodic  boundary  conditions,  because  it  obviously  does 
not  matter  where  on  the  period  the  “initial”  domain  5'(t[/^)  is  located  for 
every  j  =  0, 1,  2, ... .  Moreover,  while  the  most  universal  formulation  would 
imply  choosing  the  same  period  Z  for  all  the  coordinates,  in  the  case  when 
motion  (9)  is  characterized  by  a  predominant  direction,  the  periods  in  the 
directions  orthogonal  to  that  can  be  chosen  smaller.  For  subsequent  analysis 
in  the  current  paper,  the  periodic  setup  will  always  be  assumed. 


Next,  let  us  recast  each  formula  (21)  for  the  discrete  case  (subscript  “h”)  in 
the  form  of  a  difference: 


m2 


Ph{x,t)  =  Y. 


m2 


Uh{x,t)  =  Y 

j=mi 


j=mi 


m2  ,  mi-1  ,  and  (24) 

=  YPh\^,t)-  E  Ph\^,t)  =Y'^h\^^'^)  -  E 


j=0 


j=0 


j=0 


j=0 
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Existence  of  the  upper  limit  j  =  m2  in  the  summation  (21)  or  (24)  is  due  to 
the  causality  which  is  always  a  factor,  and  has  nothing  to  do  with  the  lacunae. 

m2  /  -N  m2  /  -N 

Therefore,  each  minuend  in  formulae  (24),  X)  Ph  oi’  Z)  could 

j=0  j=0 

have  simply  been  obtained  by  a  straightforward  time- marching  of  system  (10) 
on  the  interval  [0,  t)  in  the  foregoing  periodic  setting,  with  absolutely  no  regard 
to  either  the  partition  (16)  or  split  systems  (18).  The  full  quantity,  ph{x,t)  or 
Uh{x,  t),  cannot,  of  course,  be  obtained  by  only  marching.  To  properly  address 

mi  — 1  .  ..  mi  — 1  /  ... 

the  presence  of  the  subtrahends  Z  Ph  {x,t)  and  Z  (a;,  t)  in  formulae 

j=0  j=0 

(24),  let  us  hrst  symbolically  write  down  the  time-marching  scheme  that  would 
apply  to  system  (10),  as  well  as  to  all  systems  (18): 

Ph(x,t  +  At)  =  v(ph(x,t),Uh(x,t),q^oi(x,t)), 

(25) 

Uh(x,t  +  At)  =u(^ph(x,t),Uh(x,t),  b^oi(x,t)y 


Scheme  (25)  is  chosen  two-level  explicit  for  simplicity  only,  this  is  by  no  means 
a  limitation,  and  the  analysis  for  multi-level  schemes  can  be  found  in  [12]. 
Consider  now  a  particular  moment  of  time  t  that  corresponds  to  the  change 
in  the  lower  summation  limit  in  formulae  (21),  and  accordingly  (24),  from 
j  =  mi  to  j  =  mi  -|-  1,  i.e.,  such  t  that 


mi  -|-  1 


(26) 


Combining  formulae  (24)  and  (25),  we  will  then  have: 

Ph{x,t  +  At)  =  V(ph{x,t),  Uh{x,t),q^o\{x,t))  -p^^^\x,t  +  At), 

(27) 

Uh{x,t  +  At)  =u{ph{x,t),  Uh{x,t),  b^oi{x,t)'^  -  ul"^^\x,t  +  At). 

In  other  words,  when  the  current  moment  of  time  t  satishes  the  “switch” 
condition  (26),  the  terms  p^^^\x,t  At)  and  u]^^\x,t  +  At)  need  to  be 
explicitly  subtracted  from  the  respective  overall  expressions,  see  (27),  on  top 
of  the  standard  time-marching  step  as  per  (25).  This  basically  amounts  to 
the  required  change  of  the  upper  summation  limit  in  both  subtrahends  of 
formulae  (24)  from  mi  —  1  to  mi.  We  will  also  assume  hereafter  that  similarly 
to  the  original  differential  equations  (10),  the  scheme  (25)  will  satisfy  the 
linear  superposition  principle.  Then,  the  next  time  step  after  the  one  dehned  by 
formulae  (27),  i.e.,  the  step  t  +  At  1 — >■  t  +  2At,  shall  only  be  done  by  marching 
(25).  Indeed,  the  subtracted  quantities  p^ff'^\x,t  -\-  At)  and  ujf^^\x,t  At) 
will  carry  over  to  all  the  steps  that  follow  (27)  due  to  the  linearity.  The  genuine 
“unperturbed”  marching  can  thus  continue  till  the  next  switching  moment  t, 
i.e.,  till  condition  (26)  is  satished  by  mi  -|-  2  and  mi  -|-  1.  At  this  moment,  the 
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quantities  ,t  +  /S.t)  and  +  At)  will  need  to  be  subtracted, 

and  then  the  procedure  will  cyclically  repeat  itself. 

Thus,  the  lacunae-based  algorithm  can  be  implemented  as  a  conventional  time¬ 
marching  procedure  supplemented  by  repeated  subtraction  of  the  retarded 
terms.  The  subtraction  moments  are  known  up-front  and  separated  from  one 
another  by  equal  time  increments.  The  subtracted  terms  t  +  At)  and 

uj^^\x,t  +  At)  are  legitimately  called  retarded  because  for  a  given  moment 
of  time  t  that  satishes  (26),  they  are  generated  by  the  RHSs  and 
that  are  active  in  the  past,  on  the  time  interval  +  T]  =  [{ami  — 

\)T,  {ami  +  \)T]  =  [t  —  Tint,^  —  ^int  +  T].  Of  course,  the  actual  subtracted 
quantities  p^^^\x,t  +  At)  and  uj^^\x,t  +  At)  need  to  be  re-computed  for 
every  mi  independently  of  the  primary  time-marching  procednre.  This  is  done 
by  means  of  the  same  scheme  (25)  applied  to  the  corresponding  system  (18). 

To  conclude  this  section,  we  note  that  the  original  idea  behind  using  the 
lacunae  is  to  keep  the  number  of  terms  in  sums  (21)  or  (24)  fixed  and  non¬ 
increasing,  while  still  guaranteeing  that  the  solution  for  x  &  S{t)  will  be  the 
same  as  if  we  integrated  system  (10)  continuously  starting  from  t  =  0.  Besides, 
existence  of  the  lower  summation  limit  j  =  mi  in  (21)  and  (24),  i.e.,  repeated 
snbtraction  of  the  retarded  terms,  serves  an  additional  important  pnrpose.  It 
keeps  the  reflected  waves  from  coming  back  into  the  domain  S{t)  after  the  time 
interval  Tint  has  elapsed  since  the  inception  of  any  given  component  (16)  of 
the  RHS.  Unless  explicitly  snbtracted,  these  waves  generated  by  the  sources 
and  on  the  time  interval  +  T]  for  every  given  j,  will  start 

“contaminating”  the  solution  on  S{t)  right  after  the  moment  -|-  T^t. 


5  Lacunae-Based  ABCs 


Suppose  that  the  original  formulation  of  the  problem  that  we  want  to  solve 
involves  the  entire  space  bnt  we  are  only  interested  to  hnd  a  fragment  of 
the  overall  solntion  dehned  on  the  domain  S{t).  As  in  Section  4,  the  latter  is 
snpposed  to  have  a  hxed  shape  and  hnite  size,  but  is  allowed  to  move  according 
to  the  law  (9).  While  not  making  any  specific  assnmptions  regarding  the  natnre 
of  the  phenomena/processes  that  are  going  on  inside  S{t),  we  assnme  that 
ontside  S{t),  i.e.,  Vt  >  0  and  x  ^  S{t),  the  appropriate  model  would  be  based 
on  the  homogeneous  acoustics  system: 


1  dp{x,t) 
dt 


+  poV  ■  u{x,t) 


du{x,  t) 


-h  Vp{x,t) 


0, 

0. 


(28) 
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We  assume  that  the  overall  problem,  i.e.,  the  interior  one  that  we  do  not 
specify,  combined  with  the  exterior  one,  which  is  governed  by  system  (28),  is 
uniquely  solvable  on  In  other  words,  our  model  may  include  some  possibly 
complex  phenomena  conhned  to  the  bounded  region  S{t)  that  manifest  them¬ 
selves  by  the  acoustic  sound  in  the  far  held.  The  objective  is  to  actually  solve 
the  problem  on  the  domain  S{t)  using  a  numerical  method,  but  truncate  all 
of  its  exterior  replacing  it  with  the  ABCs  on  the  external  boundary  of  S{t). 
The  ideal  or  exact  ABCs  would  make  the  foregoing  replacement  equivalent, 
which  means  that  the  solution  obtained  this  way  would  coincide  on  S{t)  with 
the  corresponding  fragment  of  the  original  inhnite-domain  solution.  The  lat¬ 
ter,  however,  is  not  actually  available  because  the  far-held  sound  propagation 
governed  by  (28)  cannot  be  computed  directly  at  an  acceptable  cost. 

From  the  viewpoint  of  an  observer  inside  S{t),  the  role  of  the  ABCs  at  dS{t)  is 
only  to  guarantee  that  this  bonndary  will  behave  exactly  as  if  the  domain  S{t) 
were  surronnded  by  an  inhnite  linear  isotropic  sonnd-condncting  medinm.  In 
particular,  the  boundary  dS(t)  may  not  rehect,  fnlly  or  partially,  any  ontgoing 
waves.  Therefore,  as  far  as  the  ABCs  are  concerned,  we  indeed  do  not  have 
to  be  very  specihc  regarding  the  nature  of  the  problem  inside  S{t),  provided 
that  the  overall  interior /exterior  problem  does  have  a  nniqne  solntion. 

The  latter  assumption  is  of  central  importance.  However,  justifying  it  for  every 
particular  formulation  is  beyond  the  scope  of  the  cnrrent  paper.  For  example, 
the  overall  problem  may  be  linear,  and  therefore,  nniqnely  solvable,  snch  as 
acoustic  scattering  from  a  given  solid  inside  S{t).  On  the  other  hand,  one  may 
consider  a  snbstantially  more  complex  model  inside  S(t),  snch  as  the  unsteady 
flow  aronnd  a  manenvering  aircraft.  The  flow  linearizes  in  the  far  held  thus 
reducing  the  original  Euler’s  equations  to  system  (28)  at  a  distance  from  the 
aircraft.  In  this  case,  justifying  the  overall  solvability  theoretically  is  difficult 
at  best.  Besides,  linearization  in  the  far  held  is  only  an  approximation  that 
becomes  more  accurate  the  further  away  from  the  aircraft  it  is  introduced. 
This  is  going  to  ahect  the  hnal  accnracy  of  the  resulting  ABCs.  ^  Hereafter, 
we  will  not  be  primarily  concerned  with  either  the  validity  of  the  linear  model 
(28)  ontside  S(t),  or  with  the  overall  solvability  issne.  We  will  rather  focus  on 
constrncting  the  ABCs  nnder  these  assnmptions  keeping  in  mind  that  they 
may  need  to  be  corroborated  independently  for  each  specihc  case. 

There  will  be  two  stages  in  constrncting  the  ABCs.  First,  an  auxiliary  prob¬ 
lem  will  be  formulated  that  will  have  the  exact  same  solution  ontside  S{t)  as 
the  original  combined  problem  does,  bnt  will  be  linear  thronghont  the  entire 
space  and  will  be  driven  by  the  specially  derived  sonrce  terms  concentrated 
only  inside  S{t).  In  other  words,  this  problem  will  be  of  the  type  considered 


^  Our  previous  steady-state  ABCs  for  external  flows  [22]  were  based  on  the  far-held 
linearization  and  have  still  proven  superior  to  other  methods. 
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in  Section  4.  Next,  the  auxiliary  problem  will  be  solved  on  a  domain  slightly 
larger  than  S(t)  using  the  lacunae-based  methodology  of  Section  4,  and  its 
solution  obtained  right  outside  S(t)  will  be  used  to  supply  the  boundary  con¬ 
ditions  for  the  original  interior  problem  solved  inside  S{t).  In  practice,  the 
two  aforementioned  stages  will  be  meshed  together  so  that  both  the  interior 
problem  and  the  auxiliary  problem  are  time-marched  concurrently.  The  entire 
algorithm  will  be  implemented  directly  on  the  discrete  level. 

We  note  that  the  intention  to  use  lacunae  for  solving  the  auxiliary  problem 
imposes  a  certain  restriction  on  the  class  of  admissible  formulations,  because 
by  doing  so  the  acoustic  far  held  is  assumed  vorticity  free.  It  is  also  known, 
however,  that  for  the  linearized  hows  the  vortical  and  acoustic  modes  essen¬ 
tially  decouple.  This  suggests  that  the  proposed  methodology  can  potentially 
be  used  to  set  the  ABCs  for  the  acoustic  part  only,  whereas  for  vorticity  a 
diherent  method  may  be  employed  (convection  along  entropy  characteristics). 

Let  us  dehne  a  subdomain  Ss(t)  C  S(t)  such  that  x  G  S^it)  if  and  only 
if  a;  G  S(t)  and  dist{x,  dS(t))  >  e,  where  e  >  0,  and  introduce  a  multiplier 
function  fx  =  ^{x,  t)  that  is  smooth  across  the  entire  space  and  Vf  >  0  satishes: 

fO,  x^Ss{t), 

/r(a;,t)  =  Jl,  x  ^  S{t),  (29) 

iG(0,l),  X  e  S{t)\Se{t), 

The  curvilinear  strip  S(t)\Se(t)  of  width  e  adjacent  to  the  boundary  dS(t)  of 
the  domain  S{t)  from  inside  will  hereafter  be  called  the  transition  region. 

Assume  now  that  the  solution  to  the  combined  interior /exterior  problem  is 
known  on  for  f  >  0.  Clearly,  it  should  satisfy  p{x,  0)  =  0  and  ^(0;,  0)  =  0 
for  X  G  M^\S'(0).  It  is  also  important  to  realize  that  the  unknown  quantities 
inside  and  outside  S(t)  do  not  necessarily  have  to  be  the  same.  For  example,  if 
the  interior  problem  is  that  of  the  flow  around  an  aircraft,  then  the  unknowns 
inside  S{t)  will  be  the  actual  flow  variables,  whereas  the  variables  p{x,t) 
and  u{x,t)  in  system  (28)  that  is  used  outside  S(t)  are  perturbations  with 
respect  to  the  corresponding  background.  We  can,  however,  assume  with  no 
loss  of  generality  that  the  dehnitions  of  the  interior  and  exterior  quantities 
are  equivalent  on  the  (narrow)  transition  region  S{t)\Se{t).  In  the  foregoing 
example  with  an  aircraft,  this  assumption  would  imply  that  linearization  in 
the  transition  region  is  still  valid. 

Having  re-de£ned  the  interior  solution  in  terms  of  the  exterior  quantities  on 
S(t)\Ss(t),  we  multiply  the  overall  solution  everywhere  by  p{x,t)  of  (29): 

p{x,t)\ — ^  p{x,t)p{x,t)  =p{x,t), 

u{x,t)  I — >•  fi{x,t)u{x,t)  =  u{x,t). 
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Fig.  1.  Schematic  geometry  of  the  auxiliary  sources  region. 

We  emphasize  that  in  order  to  obtain  p{x,t)  and  u{x,t),  see  (30),  we  do 
not  need  to  know  p{x,t)  and  u{x,t)  on  Ss(t),  because  the  multiplier  p  is 
equal  to  zero  there  anyway.  Moreover,  multiplication  (30)  will  not  change  any 
quantities  on  MP\S{t).  With  that  in  mind,  we  apply  the  differential  operator 
from  the  left-hand  side  of  system  (28)  to  the  modihed  solution  p{x,  t),  u{x,  t), 
see  (30),  and  obtain: 


^  ^  S{t)\Se{t), 
dt  [pogvoi,  X  e  S{t)\Ss{t), 


^  X  ^  S{t)\Seit), 

''“ai  ’’  x€S{t)\s.{t). 


(31) 


The  result  will  indeed  be  zero  on  Se(t),  see  (31),  because  n{x,t)  =  0  for  a;  G 
Se(t).  The  result  is  also  zero  on  the  entire  exterior  domain  MP\S(t),  see  (31), 
because  system  (28)  holds  there,  and  the  quantities  p(a;,  t)  and  u{x,  t)  coincide 
with  p{x,t)  and  u{x,t),  respectively,  on  MP\S(t),  as  the  second  identity  (29) 
combined  with  dehnitions  (30)  suggest.  The  result  in  (31)  may  differ  from  zero 
only  in  the  transition  region  S(t)\Se(t),  where  we  actually  generate  the  RHSs 
qmi{x,t)  and  b^oi{x,t).  Note,  the  smoothness  of  the  original  solution  p{x,t), 
u{x,t),  as  well  as  that  of  the  multiplier  p{x,t),  guarantee  that  these  new 
auxiliary  sources  and  b^oi{x,t)  will  be  smooth  compactly  supported 

functions.  In  Figure  1,  we  schematically  depict  the  geometry  of  the  region  on 
which  the  auxiliary  sources  are  dehned. 


Next,  we  substitute  the  auxiliary  sources  q^oi{x,t)  and  b^oi{x,t)  of  (31)  into 
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the  RHS  of  system  (10)  and  integrate  the  latter  subject  to  the  homogeneous 
initial  conditions.  Because  of  the  overall  regularity  that  we  have  assumed  ahead 
of  time,  the  solution  of  this  new  problem  that  we  will  further  refer  to  as  the 
auxiliary  problem,  will  be  unique  and  will  coincide  with  the  modihed  functions 
p{x,t)  and  u{x,t)  of  (30)  on  the  entire  space.  This  means  that  on  MP\S{t) 
the  solution  to  the  auxiliary  problem  will  coincide  with  the  original  exterior 
solution  p{x,  t),  u{x,t).  We  emphasize  that  we  do  not  need  to  explicitly  know 
the  original  exterior  solution  on  MP\S{t)  in  order  to  obtain  the  source  terms 
qmi{x,  t)  and  6voi(®,  t)  that  drive  the  auxiliary  problem,  we  only  need  to  know 
that  it  satishes  the  homogeneous  acoustics  system  (28). 

Altogether,  we  have  split  the  original  problem  into  two:  The  linear  auxiliary 
problem  that  needs  to  be  solved  on  the  entire  space,  and  the  interior  problem 
on  S{t)  that  will  be  integrated  with  the  external  boundary  data  provided  by 
the  solution  of  the  auxiliary  problem.  However,  to  apply  the  lacunae-based 
algorithm  of  Section  4  to  the  auxiliary  problem,  it  may  need  to  be  modihed. 
Namely,  conservativeness  of  the  auxiliary  forcing  6voi  of  (31)  must  be  main¬ 
tained  so  that  to  guarantee  existence  of  the  lacunae  in  the  auxiliary  solutions, 
see  Propositions  1  and  2.  Let  us  assume  that  in  the  original  combined  problem 
the  velocity  potential  does  exist  in  the  solutions  of  system  (28)  on  MP\S(t). 
It  is  reasonable  to  think  that  it  will  exist  on  the  transition  region  S{t)\S^{t) 
as  well,  and  that  altogether  the  potential  function  (p{x,t)  will  be  smooth. 
However,  multiplication  (30)  may  ruin  the  conservativeness.  As  such,  we  can 
hrst  reconstruct  the  velocity  potential  ip{x,t)  for  x  G  S(t)\S^(t)  based  on 
the  computed  interior  quantities,  multiply  it  by  /i,  and  only  then  obtain  the 
modihed  velocity  vector  u{x,t),  thus  replacing  (30)  with  the  following  steps: 

p{x,t)  I — ^  p{x,t)p{x,t)  =p{x,t), 
u{x,t)  =  \/(p{x,t),  X  e  S{t)\S,{t), 

(p{x,t)  I - >•  fi{x,t)(p{x,t)  =  (p{x,t), 

u{x,t)  I — >  \/(p{x,t)  =  u{x,t). 

Subsequently,  the  modihed  functions  p{x,  t)  and  u{x,t)  of  (32)  are  substituted 
into  (31)  to  obtain  the  auxiliary  RHSs. 

The  version  of  the  algorithm  that  we  actually  implement  in  the  current  pa¬ 
per  does  include  reconstruction  of  the  potential  in  accordance  with  (32).  We 
realize,  of  course,  that  this  is  by  no  means  a  must.  In  fact,  all  we  need  is 
a  smooth  extension  of  the  exterior  solution  inwards  that  would  transition  to 
zero  at  a  “depth”  e,  and  would  also  maintain  conservativeness  of  the  velocity 
held.  There  may  be  diherent  ways  of  obtaining  this  extension,  not  necessarily 
based  on  applying  p  to  the  interior  solution  on  the  transition  region. 

To  set  the  discrete  ABCs  on  the  boundary  dS(t),  we  will  need  to  apply  the 
lacunae-based  algorithm  for  solving  the  auxiliary  problem  (10),  (31),  (32)  on 
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a  domain  slightly  larger  than  S(t).  We  take  5  >  0  (to  be  specihed  later)  and 
dehne  Ss{t)  =  {(a;,  t)|  dist  (a;,  S'(t))  <  5};  clearly,  S^{t)  C  S{t)  C  Ss{t),  see 
Figure  1.  We  also  replace  diamS'(t)  =  d  by  diamSsit)  =  d  +  2S  in  the  formulae 
for  the  integration  interval  and  auxiliary  domain  size  [cf.  (20),  (13)]: 


T. 


int 


{d  +  26)+T{c  +  k) 
c  —  k 


Z  —  (d  +  2d)  +  (c  +  k)Ti^t.  (33) 


Next,  assume  that  there  is  a  space-time  grid  N  x  T,  on  which  a  discrete  ap¬ 
proximation  to  the  auxiliary  problem  is  built.  The  spatial  grid  N  consists  of 
the  nodes  n  in  the  three-dimensional  space,  whereas  the  temporal  grid  T  is 
composed  of  the  time  levels  I  =  0, 1,  2, ... .  The  grid  N  is  actually  introduced 
on  the  auxiliary  domain  of  size  Z  given  by  (33),  and  periodic  boundary  con¬ 
ditions  are  imposed.  As  no  grid  adaptation  is  needed,  it  is  most  convenient  to 
simply  use  a  uniform  Cartesian  grid.  We  also  note  that  the  original  problem 
solved  inside  S{t)  does  not  have  to  be  approximated  on  the  same  grid.  In  the 
most  general  situation,  we  will  have  different  grids  for  the  interior  problem  and 
for  the  exterior /auxiliary  problem.  Then,  in  the  transition  region  S{t)\Se{t), 
where  the  dehnitions  of  the  unknown  quantities  for  both  problems  are  equiv¬ 
alent,  we  may  need  to  employ  a  chimera-type  grid  strategy,  i.e.,  interpolate 
in-between  the  overlapping  grids.  For  the  analysis  in  the  current  paper,  how¬ 
ever,  we  will  simply  assume  that  the  quantities  from  the  interior  solution  are 
already  available  on  the  grid  sub-domain  {N  x  T}  fl  {S(t)\Se(t)}. 

Let  us  denote  by  N*,  I  =  0,1,2,...,  the  corresponding  time  levels  of  the 
grid  N  X  T,  and  by  the  stencil  of  the  scheme  associated  with  the  node 
(n,  /)  G  N  X  T.  For  simplicity,  we  will  assume  that  the  auxiliary  problem  (10), 
(31),  (32)  is  time-marched  by  an  explicit  scheme,  and  that  the  node  (n,  /)  G 
is  the  actual  upper-level  node  on  the  q  +  1-level  stencil,  i.e.,  fl  {N  x  T}  C 
U  U  . . .  U  and  fl  =  (n,  /).  In  Section  4,  we  have  assumed 

g  =  1,  see  formula  (25).  Denote  by  the  sub-grid  of  that  belongs  to  the 
interior  domain,  i.e.,  fl  where  is  the  l-th  time  level  in  actual 

units  (“seconds”),  in  the  simplest  case  =  I  At.  Then,  introduce  the  sum  of 
the  interior  sub-grids  for  all  time  levels:  N+  =  U  N]/  U  U  . . .  C  N  x  T. 
Finally,  consider  a  somewhat  larger  sub-grid  of  N  x  T:  N_|_  =  U 

(n,0eN+ 

which  is  simply  a  composition  of  all  the  stencils  obtained  when  the  upper- 
level  node  (n,  /)  sweeps  the  grid  N+;  clearly,  N_|_  C  N+.  The  part  of  the  grid 
N+  that  does  not  belong  to  N+  is  called  the  grid  boundary  and  is  denoted 
7  =  N+\N+.  We  will  require  that  the  domain  Ssit’')  be  chosen  so  that  on 
every  time  level  /  =  0, 1,  2, . . . ,  all  of  the  grid  belong  to  this  domain: 
N+  C  Ss(t’‘);  equivalently,  we  may  require  that  7*  C  Ss(t’‘).  We  note  that  the 
grid  boundary  7  is  a  narrow  fringe  of  grid  nodes  that  follows  the  geometry  of 
dS{t).  Therefore,  the  parameter  S  may  be  chosen  small,  on  the  order  of  a  few 
grid  sizes  depending  on  the  specihc  structure  of  We  also  note  that  the  way 
we  have  introduced  the  grid  boundary  7  is  actually  a  simplihed  interpretation 
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of  the  rigorous  general  construction  that  is  a  part  of  the  dehnition  of  discrete 
Calderon’s  potentials  and  boundary  projection  operators,  the  latter  are  used 
in  [14]  as  a  universal  apparatus  for  setting  the  ABCs. 

We  will  now  describe  one  time  step  of  the  combined  time-marching  algorithm 
that  involves  the  lacunae-based  ABCs.  We  will  assume  that  all  time  steps  are 
identical  and  as  such,  will  provide  an  inductive  description  of  the  algorithm. 

Suppose  that  we  have  obtained  the  solution  for  up  to  a  given  time  level  1. 
This  means  that  the  solution  is  known  not  only  on  the  interior  domain,  but 
also  on  the  grid  boundary  —  at  the  levels  7^ . . .  ,7^”'^  that  are  immediately 
needed  for  advancing  the  next  time  step,  as  well  as  at  all  the  preceding  levels. 
In  particular,  one  may  think  about  starting  the  computation  from  the  known 
(homogeneous)  initial  conditions.  First,  we  make  one  interior  time  step  and 
obtain  the  discrete  solution  everywhere  inside  including  the  transition  region, 
i.e.,  the  grid  area  r]{S ,  where  the  solution  is  assumed  to  be 
dehned  in  terms  of  the  exterior  acoustic  quantities  p{x,t)  and  u{x,t).  Then, 
we  perform  the  modihcation  (32)  in  the  discrete  framework,  which  involves 
reconstruction  of  the  velocity  potential  on  the  grid  in  the  transition  region 
A  straightforward  approach  to  that  is  contour  integration 
along  the  grid  lines;  it  has  proven  quite  robust  in  our  simulations,  see  Section  6, 
taking  in  to  account  that  the  potential  only  needs  to  be  reconstructed  on  a 
narrow  near-boundary  strip.  Having  gotten  the  modified  quantities  (32)  on  the 
grid  sub-domain  n{S ,  we  apply  the  discrete  version  of  (31) 
and  obtain  one  more  time  level  of  the  discrete  auxiliary  RHSs.  If  the  scheme 
written  on  the  stencil  approximates  system  (10)  with  the  design  accuracy 
at  some  node  (n  —  no,l  —  /q)  G  (clearly,  Iq  <  q),  then  the  discrete  RHSs 
should  be  referred  to  this  same  node  (n  —  Uq,  I  —  Iq)  and  as  such,  advancing 
the  interior  solution  till  (/  -|-  1)  would  mean  building  the  auxiliary  RHSs  till 
{I  +  1  —  Iq) .  Next,  we  make  one  time  step  for  the  auxiliary  problem  with  these 
newly  updated  RHSs,  and  obtain  its  solution  on  fl  Ss(t^~^^).  Since  we 
have  chosen  <5  >  0  so  that  C  Ss(t’'~^^),  we  determine  that  the  solution  to 
the  auxiliary  problem  will,  in  particular,  be  available  on  7^+^.  This  concludes 
one  full  time  step,  because  once  we  know  the  solution  on  all  time  levels  up 
to  (/  -|-  1)  everywhere  including  the  grid  boundary,  we  can  advance  the  next 
interior  time  step,  etc. 

The  lacunae-based  methodology  for  solving  the  auxiliary  problem  includes 
cyclic  subtractions  (27)  of  the  retarded  terms  on  top  of  the  straightforward 
time-marching.  It  is  important  to  realize  that  once  a  particular  component  has 
been  subtracted,  there  will  never  be  a  need  to  analyze/incorporate  it  again 
in  the  course  of  computation.  In  other  words,  the  subtracted  terms  are  com¬ 
pletely  disregarded  from  the  moment  of  subtraction  further  on.  Therefore, 
the  corresponding  partition  elements  (16)  of  the  auxiliary  RHSs  can  be  dis¬ 
regarded  as  well.  Consequently,  even  when  integrating  over  arbitrarily  long 
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time  intervals,  we  only  need  to  keep  a  finite  amount  of  the  past  information, 
namely,  the  auxiliary  RHSs  (31)  dehned  on  the  interval  of  duration  Tint,  see 
(33),  that  immediately  precedes  the  current  moment  of  time.  This  makes  the 
extent  of  temporal  nonlocality  of  the  proposed  ABCs  fixed  and  limited. 

The  proposed  ABCs  guarantee  that  the  external  artihcial  boundary  be  com¬ 
pletely  transparent  for  all  the  outgoing  waves.  Indeed,  they  simply  allow  these 
waves  to  propagate  beyond  the  boundary  and  then  prevent  reflections  from 
re-entering  the  domain  by  eliminating  the  retarded  components  of  the  solu¬ 
tion  in  a  timely  fashion.  The  ABCs-related  computer  expenses  per  unit  time 
or  per  time  step  remain  hxed  and  non-growing,  which  is  accounted  for  by  the 
lacunae-based  integration.  For  explicit  schemes,  the  operation  count  is  pro¬ 
portional  to  the  number  of  auxiliary  grid  nodes.  The  overall  actual  cost  is,  of 
course,  higher  than  it  would  have  been  if  the  integration  was  performed  on  the 
domain  S{t)  only,  because  the  auxiliary  domain  is  larger.  The  relation  between 
the  sizes  is  given  by  (33),  and  it  also  roughly  indicates  what  the  ratio  of  the 
work  may  be.  However,  we  obviously  cannot  integrate  on  S{t)  alone,  without 
boundary  conditions.  In  this  perspective,  the  overall  performance  assessment 
must  include  the  increase  of  the  cost  that  should  be  “weighted  against”  other 
characteristics  such  as  accuracy  of  the  boundary  treatment  and  the  range  of 
problems  that  can  be  analyzed. 


6  Numerical  Demonstrations 


For  our  numerical  simulations,  we  assume  axial  symmetry  and  employ  the 
(r,  z)  cylindrical  coordinates  to  account  for  the  important  three-dimensional 
effects  using  a  two-dimensional  spatial  geometry.  Let  u  =  u{r,z,t)  and 
w  =  w(r,z,t)  be  the  radial  and  axial  components  of  the  acoustic  velocity, 
respectively,  and  let  p  =  p{r,  z,  f)  still  denote  the  acoustic  pressure.  Let  us 
also  assume  that  po  =  1-  Then,  system  (10)  becomes: 


1  dp 
dt 


1  d{r-u)  dw 

r  dr  dz 

du  dp 


dt  dr 
dw  dp 


dt  dz 


g(r,^,t), 

hfir.zfi), 

bfir,z,t). 


(34a) 


On  the  axis  r  =  0  system  (34a)  changes.  Namely,  all  the  quantities  in¬ 
volved  must  be  continuous  and  bounded.  Then,  for  the  pressure,  which  is  a 
scalar  quantity,  the  axial  symmetry  (independence  on  the  polar  angle)  im¬ 
plies:  ^  =  0.  For  the  velocity,  which  is  a  vector  quantity,  we  obtain 

=  0.  Next,  using  Taylor’s  expansion  for  r  1, 


r=0 


m(0,  zfi)  =  0  and  ^ 


r=0 
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we  have:  u{r,  •)  =  m'(0,  •)r  +  o(r)  and  consequently,  =  1  +  o(i)^ 

which  means  that  =2^  .  Therefore,  for  r  =  0  system  (34a) 

r  or  ^_Q  or  ^_q  ’  j  \  y 

transforms  into  the  system  of  two  equations: 


1  dp 
dt 


+  2 


du 

dr 

dw 


dw 

dz 

dp 

dz 


q{0,z,t), 

bz{0,z,t), 


(34b) 


and  the  radial  momentum  equation  degenerates  for  r  =  0. 


The  domain  S{t)  will  be  a  ball  of  hxed  diameter  d  centered  on  the  axis,  with 
the  given  axial  velocity  and  coordinate  of  the  center  [cf.  formulae  (9)]: 


Wo  =  Wo(t), 


Zo  =  =  f  Wo{T)dT. 

Jo 


(35) 


Obviously,  as  the  motion  (35)  is  aligned  with  the  axis  r  =  0,  it  does  not  break 
the  axial  symmetry.  We  take  the  diameter  d  =  1.8  and  the  speed  of  sound 
c  =  1;  the  functions  Woit),  and  Zoit)  of  (35)  will  be  specihed  later. 


The  auxiliary  domain  is  a  rectangle  [0,i?]  x  [— Z/2,Z/2]  of  variables  {r,z), 
with  the  actual  sizes  R  =  tt  and  Z  =  27r.  The  boundary  conditions  for  all 
variables  are  periodic  with  the  period  Z  of  (33)  in  the  2:  direction: 

{p,  u,  w}{r,  z  ±  Z,t)  =  {p,  u,  w}{r,  z,  t).  (36a) 

In  the  radial  direction,  the  boundary  conditions  cannot  be  periodic  because 
of  the  geometry/symmetry  considerations.  At  r  =  0,  there  is  no  need  for  the 
boundary  conditions  at  all;  instead,  we  have  system  (34b)  and  u{Q,zR)  =  0. 
At  r  =  R,  the  boundary  conditions  need  to  be  provided  for  the  homogeneous 
counterpart  of  system  (34a),  and  we  set: 

p{R,z,t)  =  0,  w{R,z,t)  =  0.  (36b) 

dr  r=R 

We  note  that  the  reflecting  properties  of  boundary  conditions  (36a)  and  (36b) 
are  basically  immaterial  for  reconstructing  the  inhnite-domain  solution  on 
S(t),  as  long  as  the  geometric  restrictions  discussed  in  the  beginning  of  Sec¬ 
tion  4  are  honored.  We  only  need  to  make  sure  that  the  auxiliary  problem  is 
uniquely  solvable,  which  boundary  conditions  (36a),  (36b)  do  provide  for.  We 
leave  out  a  detailed  justihcation  of  the  latter  statement,  only  mention  that  the 
homogeneous  Robin  boundary  condition  (36b)  for  u  follows  from  the  homo¬ 
geneous  Dirichlet  boundary  conditions  (36b)  for  p  and  w  combined  with  the 
homogeneous  counterpart  of  the  continuity  equation  in  system  (34a). 


The  acoustic  test  solution  that  we  will  be  validating  our  numerical  method 
on  is  supposed  to  have  a  conservative  velocity  held.  Therefore,  it  will  be  con- 
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venient  to  first  construct  the  test  solution  for  the  potential,  and  then,  by 
differentiating  it,  obtain  the  acoustic  quantities,  see  (6).  The  latter  will  sub¬ 
sequently  be  reconstructed  by  integrating  numerically  the  acoustics  system, 
and  continuous  and  discrete  quantities  will  be  compared  against  one  another. 


According  to  Proposition  1,  velocity  potential  must  satisfy  the  wave  equation. 
For  the  moment,  let  us  assume  that  this  wave  equation  is  driven  by  a  moving 
point  source  with  the  amplitude  y  =  x(t)  [cf.  equation  (la)]: 

-  Xo(t))  =  f(x,  t).  (37) 

Solution  of  equation  (37)  can  be  obtained  as  convolution  with  the  fundamental 
solution  of  the  wave  equation. 


_  H(t)  (5(|a;|  -  ct) 

Att  t  ’ 

see,  e.g.,  [15],  where  H{t)  is  the  Heaviside  function,  and  5(|a;|  —  ct)  is  a  single 
layer  of  unit  magnitude  on  the  expanding  sphere  of  radius  ct  [cf.  formula  (2)]: 


1 

(f  =  S  *  f  =—  dr 

An  Jo 


5(|a;  -^1  -c(t-r)) 


t  —  T 

1  6{\x  -  Xq{t)\  -  c(t  -  t)) 


x{r)S{$,  -  Xo{T))d$, 


An  Jo 
1  r’^ 


x{r)dT 


t  —  T 

\x  —  Xq{t)\5{i'  —  ct) 


(38) 


An  J\x\  {c\x  -  Xo{t)\  -  {x-  cco(r),  ffo(r)))(t  -  r) 
1  cx(r) 


X{r)du 


Anc\x  -  iCo(r)|  -  {x  -  Xo{t),  Uo{t)) 


l'=Ct 


In  formula  (38),  (• ,  ■)  denotes  the  dot  product.  For  the  integration,  we  have 
introduced  a  new  variable  v  =  \x  —  Xq{t)  \  +  cr  as  done  in  [23,  Chapter  7]. 
Evaluating  the  last  expression  in  (38)  for  v  =  ct  requires  solving  equation 


X  —  iCo(T)|  -|-  cr  =  cf 


(39) 


with  respect  to  r.  Solution  of  equation  (39)  determines  the  retarded  moment  of 
time,  at  which  the  trajectory  of  the  point  source  intersects  the  lower  portion  of 
the  characteristic  cone  with  the  vertex  {x,t).  For  the  case  of  a  straightforward 
uniform  motion,  equation  (39)  is  quadratic,  and  its  solution  r,  once  substituted 
in  (38),  yields  (p{x,t)  that  can  also  be  obtained  via  the  Lorentz  transform,  as 
done  in  [12,13].  In  general  however,  one  should  not  expect  to  be  able  to  solve 
the  nonlinear  equation  (39)  in  the  closed  form. 

Let  us  now  take  into  account  the  cylindrical  symmetry  and  straightforward 


23 


motion  (35).  Then,  formulae  (38)  and  (39)  reduce  to 


J- _ cx{t) 

-  {z 


zo{t))wo{t) 


v=ct 


(40) 


and 


^Jr'^  +  {z  -  Zq{t))‘^  +  ct  =  ct, 


(41) 


respectively.  Assuming  that  ki  <  Wo(t)  <  k2,  where  ki  and  k2  are  known,  we 
solve  equation  (41)  by  Newton’s  method  with  the  initial  guess  Tq  given  by  the 
solution  of  the  quadratic  equation  r‘^  +  {z  —  ^{ki  +  k2)ToY  =  c(t  —  tqY  that 
satishes  tq  <  t.  The  actual  law  of  motion  that  we  specify  is  [cf.  formulae  (35)]: 


Zo{t) 


3 

+ 


where  [  ■  ]  denotes  the  integer  part,  as  before,  and  {  ■  }  denotes  the  fractional 
part.  The  motion  (42)  basically  consists  of  repeated  acceleration/deceleration 
cycles  of  duration  10,  so  that  during  each  cycle  the  source  travels  a  total 
distance  of  1  along  the  2:  axis.  Both  the  velocity  Wo(t)  =  ZQ(t)  and  the  acceler¬ 
ation  ao(t)  =  z'^it)  determined  by  (42)  are  continuous  functions  of  time,  and 
0  =  ki  <  wo(t)  <  k2  =  0.1875,  i.e.,  the  subsonic  condition  is  met. 


Solution  (40)  for  the  potential  is  singular,  and  cannot  be  used  directly  to  derive 
the  acoustic  quantities  needed  for  testing  the  numerical  procedure.  To  remove 
the  singularity,  we  introduce  a  new  smooth  function  G  =  G{r)  of  the  variable 
r  =  f(r,z,t)  =  +  (z  —  Zo(t)y  such  that  G{r)  =  1  for  r  >  ^,  where 

K  <  1,  and  also  such  that  at  least  G(0)  =  G'(0)  =  G"(0)  =  G'"(0)  =  0.  Then, 
the  new  function  ip{r,z,t)  ■  G{r)  is  continuous  and  bounded  everywhere.  By 
differentiating  it,  we  dehne  the  components  of  the  reference  acoustic  solution 
[cf.  formulae  (6)]: 

d 

p{r,  z,t)  =  -  —  ((^(r,  t)  ■  G(r)^ , 

u(GZ,t)  =^((p(r,z,t)  ■  G(k)),  (43) 

d 

w(r,  z,  t)  =—(ip{r,  z,  t)  ■  G(r)J , 

which  are  regular  functions  as  well.  Note  that  99  of  (40)  is  a  retarded  potential; 
therefore,  differentiation  (43)  involves  implicit  differentiation  of  r  via  (41).  At 
the  same  time,  the  dehnition  of  G(r)  does  not  involve  any  retardation. 


The  variables  p,  u,  and  w  of  (43)  satisfy  the  homogeneous  version  of  the  acous¬ 
tics  system  (34a)  for  r  >  ^,  where  we  have  taken  k  =  0.8.  Inside  the  smaller 
ball,  i.e.,  for  r  <  substitution  of  p,  u,  and  w  of  (43)  into  the  left-hand 
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side  of  system  (34a)  will  produce  the  corresponding  source  terms  q,  br-,  and 
hz-  Note,  when  removing  the  singularity  of  the  solution  we  have  required  that 
sufficiently  many  derivatives  of  G  be  equal  to  zero  at  r  =  0;  consequently, 
the  RHSs  due  to  the  quantities  (43)  substituted  into  (34a)  will  have  no  sin¬ 
gularities  either.  Altogether,  we  have  therefore  obtained  a  reference  acoustic 
solution,  which  is  regular  everywhere,  and  which  can  he  said  to  he  generated 
by  the  sources  concentrated  inside  the  moving  domain  S(t).  Outside  S(t),  our 
reference  solution  is  basically  a  system  of  unsteady  acoustic  waves  radiated  by 
and  propagating  away  from  a  moving  source.  We  are  going  to  reconstruct  this 
solution  numerically  on  the  domain  S{t),  and  set  the  discrete  lacunae-based 
ABCs  on  its  outer  boundary  dS(t)  according  to  the  methodology  of  Section  5. 


The  acoustics  system  (34a)  is  approximated  numerically  on  the  Cartesian  grid: 

Vi  =  iAr,  Ar  =  R/Nr,  i  =  0, 1, . . .  ,Nr, 

Zj=jAz,  Az  =  Z/2Nz,  j  =  0,±1,...  ,±Nz, 

using  a  second-order  staggered  hnite-difference  scheme: 
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(45a) 


Scheme  (45a)  can  be  written  for  i  >  0,  i.e.,  away  from  the  axis  of  symmetry. 

;+i 

For  i  =  0,  we  set  =  0  and  approximate  system  (34b),  which  yields: 


1  Pi  .  -  Pi  2  z+i  /+i  -  Li  i+i 

c2  At  Ar  Az 
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(45b) 


Scheme  (45a)  is  very  similar  to  the  well-known  Yee  scheme  that  was  originally 
introduced  for  solving  the  Maxwell  equations,  see  [24] .  Periodic  boundary  con¬ 
ditions  (36a)  are  re-written  on  the  grid  as  follows: 


'+5 

=  W 


l+k 


w.A  .rm  i  (46a) 


i+ANz  +  h 
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whereas  boundary  conditions  (36b)  are  discretized  as: 


i+l 


*+2 


=  0, 


Ar 


=  0, 


w 


'+1 

Nr- 


ii+5 


=  0. 


(46b) 


We  have  used  three  successively  more  hne  square  cell  grids,  Ar  =  Az,  x 
2Nz  =  64  X  128,  128  x  256,  and  256  x  512.  The  Courant  stability  constraint 
was  applied  when  selecting  the  time  step  for  the  explicit  scheme  (45a),  (45b). 
The  grid  boundary  7  was  built  as  outlined  in  Section  6,  taking  into  account 
that  scheme  (45a),  (45b)  is  staggered.  The  latter  implies  that  we  actually  have 
three  different  stencils  for  updating  the  pressure  and  two  velocity  components 
that  are  shifted  with  respect  to  one  another.  Of  course,  the  grid  boundary  is 
constructed  concurrently  with  the  actual  time-marching.  The  parameter  S  that 
is  needed  to  accommodate  the  width  of  the  grid  boundary  7  was  taken  S  = 
|Ar.  The  functions  0(t),  and  G{r)  introduced  in  order  to  guarantee 

smoothness  at  all  the  stages  of  the  derivation,  are  constructed  as  piece-wise 
polynomials  with  four  continuous  derivatives  everywhere.  In  so  doing,  the 
multiplier  see  (29),  is  also  built  as  a  function  of  r  only.  The  varying 

amplitude  y,  see  (40),  was  chosen  in  the  form  of  a  harmonic  oscillation  with 
the  frequency  three  times  that  of  the  motion  cycles  (42).  The  width  e  of 
the  transition  region  S(t)\Ss(t),  see  Figure  1,  varied  to  demonstrate  different 
aspects  of  the  algorithm  performance.  The  parameter  a  of  (14),  (15)  was 
chosen  n  =  |.  The  actual  temporal  thickness  T  of  each  partition  element 
was  calculated  “backwards”  from  (33),  considering  that  the  domain  size,  the 
maximum  motion  speed,  see  (42),  and  the  period  Z  are  known.  When  marching 
the  auxiliary  problem,  retarded  components  of  the  solution  are  subtracted 
according  to  (27).  Each  subtracted  component  is  recomputed  as  solution  to  the 
corresponding  problem  (18).  For  a  given  j,  this  problem  is  inhomogeneous  on 
the  interval  4'^^]  =  [(aj  — 1)T,  {(7j+^)T],  and  then  it  remains  homogeneous 

on  the  interval  =  [{(jj  +  |)T,  {aj  —  1)T  -|-  Tint]-  Therefore,  we  first 

explicitly  time-march  this  system  on  its  interval  of  inhomogeneity.  Then,  we 
do  the  FFT  of  the  solution  in  the  direction  and  expansion  with  respect  to  the 
corresponding  eigenfunctions  (evaluated  numerically)  in  the  r  direction,  which 
allows  us  to  advance  the  homogeneous  solution  further  till  by  simply  raising 
the  resulting  amplihcations  factors  to  the  corresponding  powers.  We  note  that 
in  so  doing  scheme  (45a)  effectively  gets  split  into  three  discrete  wave  equations 
for  the  respective  unknown  quantities,  because  in  the  cylindrical  geometry 
different  transforms  along  r  are  needed  for  different  variables.  This  is  the  only 
instance  when  reduction  to  a  set  of  independent  wave  equations  is  employed; 
and  it  is  only  necessitated  by  a  particular  choice  of  the  coordinate  system.  If  a 
genuine  three-dimensional  computation  were  conducted  on  a  Cartesian  grid, 
the  unsplit  system  could  have  been  used  all  along. 

In  Figure  2  we  present  error  for  the  acoustic  pressure  as  it  depends  on  time  on 


26 


Lacunae-based  ABCs  for  unsteady  acoustics 


(D 

> 

O) 

o 


Grid  convergence  for  the  pressure;  accelerated  source  motion 

-0.5 
-1 
-1.5 
-2 
-2.5 
-3 
-3.5 
-4 
-4.5 
-5 
-5.5 
-6 
-6.5 
-7 
-7.5 
-8 
-8.5 
-9 
-9  5 

■  0  20  40  60  80  100 

Dimensionless  time 


Fig.  2.  Grid  convergence  study  with  lacunae-based  ABCs,  e  =  lOAr. 

all  three  grids  that  we  have  employed.  The  total  integration  time  was  100^, 
i.e.,  one  hundred  times  the  interval  required  for  the  waves  to  cross  the  domain. 
The  error  was  evaluated  in  the  maximum  norm  on  the  domain  S{t).  We  see 
that  the  algorithm  provides  for  no  long-term  error  buildup,  and  also  that  it 
displays  the  design  second-order  grid  convergence.  The  oscillations  in  the  error 
prohles  on  Figure  2  have  the  period  of  10  time  units,  and  apparently  follow 
the  acceleration/deceleration  cycles  of  the  source  motion.  Error  prohles  for 
the  velocity  components  u  and  w  look  similarly  and  we  do  not  present  them. 


We  should  mention  that  the  error  curves  on  Figure  2  were  obtained  for  a  rela¬ 
tively  wide  transition  region  S{t)\S^{t),  on  the  order  of  ten  grid  cell  sizes.  Even 
though  the  actual  width  of  this  region  decreases  with  the  rehnement  of  the 
grid,  the  issue  of  how  the  width  e  affects  the  algorithm  performance  deserves 
to  be  thoroughly  addressed.  From  the  numerical  standpoint,  the  width  of  the 
transition  region  determines  how  well  the  smooth  function  n{x,t)  is  resolved 
on  the  grid,  and  as  such,  how  smooth  the  auxiliary  RHSs  will  effectively  be. 
The  latter,  in  turn,  affect  the  quality  of  the  discrete  lacunae,  i.e.,  how  sharp 
the  aft  fronts  of  the  waves  really  are  in  the  discrete  framework.  This  is  im¬ 
portant  because  every  time  a  retarded  component  is  subtracted,  see  (27),  we 
assume  that  what  is  being  subtracted  inside  the  lacuna,  i.e.,  on  the  domain 
S{t),  is  zero,  or  more  precisely,  a  small  quantity  that  converges  to  zero  with 
the  rehnement  of  the  grid;  the  convergence,  however,  hinges  on  the  smoothness 
of  the  source  terms.  Besides  having  a  potential  effect  on  the  error  behavior. 
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the  width  of  the  transition  region  also  determines  on  how  many  grid  nodes  the 
auxiliary  RHSs  are  supported.  We  remind  that  those  RHSs  basically  control 
the  extent  of  temporal  nonlocality  of  the  lacunae-based  ABCs.  The  algorithm 
requires  keeping  them  on  the  interval  of  length  Rnt,  and  as  such,  the  more 
narrow  the  transition  region  is,  the  less  additional  storage  is  needed. 

In  Figures  3,  4,  5,  and  6  we  are  showing  similar  pressure  error  prohles  for  the 
width  of  the  transition  region  e  =  8,  6,  4,  and  2  grid  cell  sizes,  respectively. 

We  observe  that  with  the  decrease  of  e  the  error  behavior  deteriorates,  which 
is  natural  to  expect.  We  also  notice,  though,  that  the  deterioration  is  more 
visible  on  the  coarser  grids,  whereas  on  the  hnest  256  x  512  grid  it  is  much 
slower.  The  256  x  512  error  prohle  is  still  practically  flat  for  e  =  8,  see  Figure  3, 
and  it  loses  only  one  binary  order  of  magnitude  for  e  =  6  over  the  entire  long 
interval  of  integration,  see  Figure  4.  Even  for  a  rather  narrow  transition  region, 
e  =  4,  see  Figure  5,  the  hnest  grid  error  only  grows  by  less  than  a  factor  of 
2  over  the  hrst  half  of  the  integration  interval,  which  is  still  hfty  times  the 
time  needed  for  the  waves  to  cross  the  domain.  The  rates  of  error  increase 
become  practically  equal  on  all  grids  only  for  the  narrowest  transition  region 
that  we  have  tried:  e  =  2,  see  Figure  6.  Having  only  two  grid  cells  in  the 
transition  region  basically  implies  that  there  is  no  smoothing  at  all,  rather 
a  sharp  truncation,  and  instead  of  n{x,t)  we  are  using  an  equivalent  of  the 
Heaviside  function  on  the  grid.  But  even  in  this  case  we  can  see  that  the  initial 
jump  of  the  error  is  much  smaller  on  the  hne  grid  than  on  the  coarser  grids. 

As  of  yet,  we  cannot  offer  a  rigorous  theoretical  explanation  of  why  the  algo¬ 
rithm  appears  more  sensitive  to  the  quality  of  the  discrete  lacunae  on  coarser 
grids  than  on  the  fine  grid.  We  can  only  qualitatively  suggest  that  it  has  to 
do  with  the  actual  magnitude  of  those  discrete  “tails”  behind  the  aft  fronts 
of  the  waves  that  are  due  to  the  “imperfections”  in  the  auxiliary  sources,  and 
that  apparently  are  still  smaller  on  hne  grids.  Altogether,  this  phenomenon 
is  certainly  benehcial,  because  hne  grids  are  needed  for  high  overall  accuracy 
anyway,  and  at  the  same  time  they  will  allow  to  maintain  high  accuracy  of 
the  boundary  treatment  for  longer  periods  of  time.  We  should  also  note  that 
for  many  practical  applications  the  actual  integration  times  will  likely  be  not 
as  long  as  those  that  we  have  used  for  the  current  proof-of-concept.  This  will 
leave  even  less  room  for  the  error  buildup  due  to  the  boundary  treatment. 


7  Conclusions 


We  have  constructed  and  tested  the  algorithm  for  setting  highly-accurate 
global  artihcial  boundary  conditions  for  the  computation  of  time-dependent 
acoustic  waves.  This  work  is  an  extension  of  our  previous  approach  that  applied 
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3.  Grid  convergence  study  with  lacunae-based  ABCs,  e  =  8Ar. 
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Fig.  4.  Grid  convergence  study  with  lacunae-based  ABCs,  e  =  6Ar. 
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Fig.  5.  Grid  convergence  study  with  lacunae-based  ABCs,  e  =  4Ar. 
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Fig.  6.  Grid  convergence  study  with  lacunae-based  ABCs,  e  =  2Ar. 
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to  the  scalar  wave  equation  [13] .  The  algorithm  is  based  on  the  presence  of  lacu¬ 
nae  (aft  fronts  of  the  waves)  in  the  three-dimensional  wave-type  solutions.  The 
ABCs  are  obtained  directly  for  the  discrete  formulation  of  the  problem  and 
can  complement  any  consistent  and  stable  hnite-difference  scheme.  In  so  doing, 
neither  a  rational  approximation  of  non-reflecting  kernels,  nor  discretization  of 
the  continuous  boundary  conditions  is  required.  The  extent  of  temporal  non- 
locality  of  the  new  ABCs  appears  hxed  and  limited,  and  this  is  not  a  result  of 
any  approximation  but  rather  a  direct  consequence  of  the  fundamental  prop¬ 
erties  of  the  solution.  The  proposed  ABCs  can  handle  artihcial  boundaries  of 
irregular  shape  on  regular  grids  with  no  htting/adaptation  needed.  Besides, 
they  possess  a  unique  capability  of  being  able  to  handle  boundaries  of  moving 
computational  domains,  including  the  case  of  accelerated  motion. 

For  the  case  of  an  acoustic  source  engaged  in  an  accelerated  motion,  we  have 
conducted  numerical  experiments  that  corroborate  the  theoretical  design  prop¬ 
erties  of  the  algorithm.  We  are  currently  unaware  of  any  other  acoustic  ABCs’ 
algorithm  in  the  literature  with  the  capability  of  handling  the  accelerated 
motion.  We  have  also  shown  experimentally  that  when  a  key  parameter  that 
characterizes  the  algorithm  changes  so  that  to  reduce  the  overall  memory  re¬ 
quirements,  the  performance  of  the  ABCs  snffers  in  an  expected  way.  However, 
the  deterioration  of  the  long-term  performance  on  fine  grids  is  much  slower 
than  that  on  coarse  grids. 

Finally,  we  should  mention  that  even  thongh  the  full  description  of  the  algo¬ 
rithm  provided  in  the  paper  does  address  a  nnmber  of  technical  issues,  its  key 
idea  is  most  straightforward.  In  one  sentence  it  can  be  formulated  as  follows: 
One  should  not  continne  the  compntation  inside  the  lacnna  once  the  solu¬ 
tion  has  become  zero  there  due  to  “natural  causes.”  Otherwise  the  error  may 
unwarrantably  build  up  on  one  hand,  and  on  the  other  hand,  the  extent  of 
the  reqnired  temporal  pre-history  of  the  solution  may  grow  un-justi£ably  high. 
The  technical  issnes  that  we  have  discnssed  relate  primarily  to  what  to  do  with 
the  waves  outside,  and  not  inside,  the  domain  of  interest.  In  the  framework  of 
the  previous  analysis,  we  simply  allow  them  to  propagate  a  certain  distance 
away  till  they  get  reflected,  and  then  set  up  the  auxiliary  domain  so  that  the 
reflected  waves  do  not  reach  the  domain  of  interest  before  it  completely  falls 
inside  the  lacnna.  This,  however,  is  by  no  means  the  only  possible  option.  In 
fact,  any  treatment  of  the  ont going  waves  that  would  prevent  the  reflections 
from  re-entering  the  domain  of  interest  before  it  falls  into  the  laeuna  will  be 
appropriate.  At  the  same  time,  introducing  an  alternative  to  the  approach 
described  in  the  current  paper  may  be  beneficial  from  the  standpoint  of  the 
overall  computational  cost.  For  example,  a  treatment  of  the  sponge  layer  type 
that  slows  down  the  outgoing  waves,  see,  e.g.,  [25,26],  may  allow  to  reduce  the 
size  of  the  auxiliary  domain.  Alternatively,  the  lacnnae-based  approach  can 
be  combined  with  a  PML-based  treatment  for  the  waves  ontside  the  domain 
of  interest,  see,  e.g.,  the  survey  papers  [27,28].  In  any  event,  linearity  has 
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to  be  maintained,  otherwise  it  will  not  be  possible  to  partition  the  problem 
similarly  to  (16),  (18),  (19).  These  combined  approaches  will  be  subject  of  a 
future  study. 


Acknowledgements 


The  author  is  very  thankful  to  the  reviewers  of  the  paper  for  their  helpful 
comments. 


References 


[1]  D.  Givoli,  Non-reflecting  boundary  conditions,  J.  Comput.  Phys.  94  (1991)  1-- 
29. 

[2]  T.  Hagstrom,  Radiation  boundary  conditions  for  the  numerical  simulation  of 
waves,  in:  A.  Iserlis  (Ed.),  Acta  Numerica,  Vol.  8,  Cambridge  University  Press, 
Cambridge,  1999,  pp.  47-106. 

[3]  S.  V.  Tsynkov,  Numerical  solution  of  problems  on  unbounded  domains.  A 
review,  Appl.  Numer.  Math.  27  (1998)  465-532. 

[4]  M.  J.  Grote,  J.  B.  Keller,  On  nonreflecting  boundary  conditions,  J.  Comput. 
Phys.  122  (1995)  231-243. 

[5]  M.  J.  Grote,  J.  B.  Keller,  Exact  nonreflecting  boundary  conditions  for  the  time- 
dependent  wave  equation,  SIAM  J.  Appl.  Math.  55  (1995)  280-297. 

[6]  M.  J.  Grote,  J.  B.  Keller,  Nonreflecting  boundary  conditions  for  time-dependent 
scattering,  J.  Comput.  Phys.  127  (1996)  52-65. 

[7]  I.  L.  Sofronov,  Artificial  boundary  conditions  of  absolute  transparency  for  two- 
and  three-dimensional  external  time-dependent  scattering  problems,  European 
J.  Appl.  Math.  9  (1998)  561-588. 

[8]  I.  L.  Sofronov,  Non-reflecting  inflow  and  outflow  in  a  wind  tunnel  for  transonic 
time-accurate  simulations,  J.  Math.  Anal.  Appl.  221  (1998)  92-115. 

[9]  B.  Alpert,  L.  Greengard,  T.  Hagstrom,  Rapid  evaluation  of  nonreflecting 
boundary  kernels  for  time-domain  wave  propagation,  SIAM  J.  Numer.  Anal. 
37  (2000)  1138-1164. 

[10]  B.  Alpert,  L.  Greengard,  T.  Hagstrom,  Nonreflecting  boundary  conditions  for 
the  time-dependent  wave  equation,  J.  Comput.  Phys.  180  (1)  (2002)  270-296. 

[11]  E.  Michielssen,  A.  Ergin,  B.  Shanker,  D.  Weile,  The  multilevel  plane  wave  time 
domain  algorithm  and  its  applications  to  the  rapid  solution  of  electromagnetic 


32 


scattering  problems:  A  review,  in:  Mathematical  and  numerical  aspects  of  wave 
propagation  (Santiago  de  Compostela,  2000),  SIAM,  Philadelphia,  PA,  2000, 
pp.  24-33. 

[12]  V.  S.  Ryaben’kii,  S.  V.  Tsynkov,  V.  1.  Turchaninov,  Long-time  numerical 
computation  of  wave-type  solutions  driven  by  moving  sources,  Appl.  Numer. 
Math.  38  (2001)  187-222. 

[13]  V.  S.  Ryaben’kii,  S.  V.  Tsynkov,  V.  I.  Turchaninov,  Global  discrete  artificial 
boundary  conditions  for  time-dependent  wave  propagation,  J.  Comput.  Phys. 
174  (2)  (2001)  712-758. 

[14]  V.  S.  Ryaben’kii,  Nonreflecting  time-dependent  boundary  conditions  on 
artificial  boundaries  of  varying  location  and  shape,  Appl.  Numer.  Math.  33 
(2000)  481-492. 

[15]  V.  S.  Vladimirov,  Equations  of  Mathematical  Physics,  Dekker,  New- York,  1971. 

[16]  I.  Petrowsky,  On  the  diffusion  of  waves  and  the  lacunas  for  hyperbolic  equations, 
Matematicheskii  Sbornik  (Recueil  Mathematique)  17  (59)  (3)  (1945)  289-370. 

[17]  R.  Courant,  D.  Hilbert,  Methods  of  Mathematical  Physics.  Volume  II,  Wiley, 
New  York,  1962. 

[18]  L.  D.  Landau,  E.  M.  Lifshitz,  Eluid  Mechanics,  Pergamon  Press,  Oxford,  1986. 

[19]  C.  L.  Morfey,  Dictionary  of  Acoustics,  Academic  Press,  San  Diego,  2001. 

[20]  L.  Ting,  M.  J.  Miksis,  Exact  boundary  conditions  for  scattering  problems,  J. 
Acoust.  Soc.  Am.  80  (1986)  1825-1827. 

[21]  D.  Givoli,  D.  Gohen,  Nonreflecting  boundary  conditions  based  on  Kirchhoff-type 
formulae,  J.  Gomput.  Phys.  117  (1995)  102-113. 

[22]  S.  V.  Tsynkov,  External  boundary  conditions  for  three-dimensional  problems 
of  computational  aerodynamics,  SIAM  J.  Sci.  Gomp.  21  (1999)  166-206. 

[23]  P.  M.  Morse,  H.  Eeshbach,  Methods  of  Theoretical  Physics.  Parts  I  and  II, 
McGraw-Hill,  Boston,  1953. 

[24]  K.  S.  Yee,  Numerical  solution  of  initial  boundary  value  problem  involving 
Maxwell’s  equations  in  isotropic  media,  IEEE  Trans.  Antennas  Propagat.  14 
(1966)  302-307. 

[25]  M.  Israeli,  S.  Orszag,  Approximation  of  radiation  boundary  conditions,  J. 
Gomput.  Phys  41  (1981)  115-135. 

[26]  S.  Kami,  Ear-field  filtering  operators  for  suppression  of  reflections  from  artificial 
boundaries,  SIAM  J.  Numer.  Anal.  33  (1996)  1014-1047. 

[27]  S.  Abarbanel,  D.  Gottlieb,  On  the  construction  and  analysis  of  absorbing  layers 
in  GEM,  Appl.  Numer.  Math.  27  (1998)  331-340. 

[28]  E.  Turkel,  A.  Yefet,  Absorbing  PML  boundary  layers  for  wave-like  equations, 
Appl.  Numer.  Math.  27  (1998)  533-557. 


33 


